--- title: Element distribution maps keywords: fastai sidebar: home_sidebar summary: "Slicing the cube " description: "Slicing the cube " nb_path: "notebooks/60_peakmaps.ipynb" ---
In order to ultimately create the much appreciated element distribution maps, we need to slice the spectral image cube at all peak locations. This can be done using the function get_cube_slices(). This function computes Gaussian peak shapes and corresponding slice indexes for each hotmax pixel in the spectral image cube. The result can be inspected with plot_cube_slices() function.
from maxrf4u import get_cube_slices, plot_cube_slices, get_peakmaps
slices, y_gauss_list = get_cube_slices('RP-T-1898-A-3689.datastack')
ax = plot_cube_slices('RP-T-1898-A-3689.datastack');
ax.set_title('Gaussian peak profiles')
plt.tight_layout()
If life were simple we could now integrate the intensities for a given peak slice to get the corresponding element map, and be done with our work! This is for example the case for our well-known iron $K_{\alpha}$ peak located at 6.40 keV colored blue in the plot above. Let's check the precise limits of energy band that we want to integrate.
ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies')
si, sj, sk = slices[10] # These are the Fe_Ka slice indices (left, mid, right) index
print(f'The Fe_Ka peak located at {x_keVs[sj]:.02f} keV ranges from {x_keVs[si]:.02f} keV to {x_keVs[sk]:.02f} keV.')
Integration into an element map is simply reading the peak slice from the data cube, summing along the spectral axis and plotting the result.
# read Fe_Ka slice
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
FeKa_slice = cube[:,:,si:sk+1].compute() # takes a few seconds...
# integrate
FeKa_map = FeKa_slice.sum(axis=2)
# and plot...
fig, ax = plt.subplots(figsize=[8, 8])
ax.imshow(FeKa_map, vmax=75); # need to clip intensity due iron speckles
ax.set_title('Simple peak slice map for iron Fe_Ka');
For now we just ignore errors due overlapping bands and take a preliminary look at all peak maps by simply integrating all peak slices use the get_peakmaps() function.
peak_maps, keV_maps = get_peakmaps('RP-T-1898-A-3689.datastack', slices, algorithm='simple')
Although it is clear that these peak maps are not perfect yet, let's take a look. Starting point for the analysis of the maps is the summary of the peak pattern puzzle from the previous section.
from maxrf4u import plot_patterns, get_patterns
ptrn_list = get_patterns(['S', 'Ca', 'K', 'Cl', 'Fe', 'Mn', 'Cu', 'Zn', 'Pb', 'Ti', 'Ba'])
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 5])
plot_patterns(ptrn_list, ax=ax)
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax, clip_vline=False)
plot_cube_slices('RP-T-1898-A-3689.datastack', ax=ax1);
ax1.plot(x_keVs, 8*y_sum, color=[0.3, 1, 0.3], label='sum spectrum')
ax1.legend();
ax.set_title('Peak pattern summary')
plt.tight_layout()
Here is the overview of all 24 peak maps.
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
titles = [f'[{i}] {x_keVs[peak_idx]:.02f} keV' for i, peak_idx in enumerate(peak_idxs)]
multi_plot(*peak_maps_histeq, titles=titles);
Unfortunately for most other peaks then iron XRF life is not so simple. There are a two complications that we need to take into account. Especially for MA-XRF scans of drawings on paper we need to deal with many spectra for locations in contain only light elements. As a consequence a significant portion of their signal is due to inelastic (Compton) and elastic scattering by hydrogen atoms in the paper. We need to correct for this Compton ridge baseline when integrating these peaks. Our second problem is the occurrence of overlapping peaks. Let's first solve the baseline estimation problem...
The technical challenge here is that the baseline which is due to to scattering of the paper substrate will vary according to the local paper thickness, but also due to the local presence of heavier elements such as lead. A layer of lead atoms on top of a paper absorbs the incoming x-rays that therefore will not reach the paper and hence not be scattered to the detector. Spectra that contain lead peaks typically have a lower scattering baseline.
This behavior is exemplified if we compare two different hotmax spectra
The slowly varying ridge with Poisson noise and element specific emission peaks can clearly be seen in both spectra.
gap_slice, y_cont = get_continuum(x_keVs, y_sum, slices, plot_this=True)
Now I need to inspect how well this continuum baseline fits different spectra.
hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
n = 4
y_hot = hotmax_spectra[n]
y_snip = bsl.smooth.snip(y_hot, 50, )[0]
gap_slice, y_cont = get_continuum(x_keVs, y_sum, slices)
gap_average = get_gap_average(y_hot, gap_slice)
y_cont_baseline = y_cont * gap_average
fig, ax = plt.subplots(figsize=[9, 4])
hma.plot_spectrum(n, ax)
ax.plot(x_keVs, y_cont_baseline, color=[0.3, 1, 0.3])
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)
ax.set_xlim([-1, 24])
ax.set_ylim([-1, 20]);
#ax.plot(x_keVs, np.log(20*mu_Pb), color='k')
ax.plot(x_keVs, 0.6*y_maxsnip, color='orange')
import xraydb
x_eVs = 1000 * x_keVs
mu_Pb = xraydb.mu_elam('Pb', x_eVs)
import pybaselines as bsl
y_maxspectrum
y_maxsnip = bsl.smooth.snip(y_maxspectrum, 30, decreasing=True)[0]
y_maxsnip_noiseline = y_maxsnip + 0.1*y_maxsnip.max()
is_gap = y_maxsnip_noiseline > y_maxspectrum
fig, ax = plt.subplots(figsize=[9, 4])
ax.plot(x_keVs, y_maxspectrum)
ax.plot(x_keVs, y_maxsnip_noiseline)
ax.fill_between(x_keVs, y_maxsnip_noiseline, where=is_gap, color='orange', alpha=0.5)
ax.fill_between(x_keVs, y_maxsnip_noiseline, where=(labeled==19), color='r', alpha=0.5)
ax.set_xlim([-1, 25])
ax.set_ylim([-2, 50]);
ax.scatter(x_keVs[hill_top_idx], y_maxsnip[hill_top_idx])
Need to inspect the width of each gap band to select best one. Or instead test dask boolean indexing.
is_test = np.zeros_like(x_keVs).astype(bool)
hill_top_idx = np.argmax(y_maxsnip)
is_test[1000:1100] = True
test_slice = cube[:,:,is_test]
test_slice
import skimage.morphology as skm
sizes[19-1]
labeled = skm.label(is_gap)
sizes = []
for l in range(1, labeled.max()):
d = np.sum(labeled == l)
sizes.append(d)
sizes
y_sumsnip = bsl.smooth.snip(y_sum, 30, decreasing=True)[0]
fig, ax = plt.subplots(figsize=[9, 4])
ax.plot(x_keVs, y_sum)
ax.plot(x_keVs, y_sumsnip)
slices, y_gauss_list = get_cube_slices('RP-T-1898-A-3689.datastack')
import matplotlib.cm as cm
y_continuum = smooth_compton(y_sum)
# peak free region for continuum normalization
# not sure how to automate this choice
# 18 - 19 is near continuum top and fairly large region
left_peak_idx = 18
right_peak_idx = left_peak_idx + 1
gap_band_i = slices[left_peak_idx][2] + 1
gap_band_j = slices[right_peak_idx][0] -1
continuum_gap_mask = np.zeros_like(x)
continuum_gap_mask[gap_band_i:gap_band_j+1] = 1
print(f'Gap region width: {gap_band_j - gap_band_i}')
y_
x = x_keVs
y_maxspectrum = ds.read('maxrf_maxspectrum')
#def plot_cube_slices(slices, x, y):
mask_list = []
for [si, sj, sk] in slices:
mask = np.zeros_like(x)
mask[si:sk+1] = 1
mask_list.append(mask)
sum_mask = np.array(mask_list).sum(axis=0)
fig, ax = plt.subplots(figsize=[9, 4])
for i, mask in enumerate(mask_list):
ax.fill_between(x, 100*mask, where=mask, alpha=0.2)
sj = slices[i][1]
#ax.axvline(x[sj], color='r')
ax.plot(x, sum_mask)
ax.plot(x, y_maxspectrum)
ax.plot(x, y_compton, label='Compton')
ax.legend()
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)
ax.fill_between(x, y_maxspectrum, where=continuum_gap_mask, color='r')
ax.set_xlim([-1, 25]);
In this case, the best
datastack_file = 'RP-T-1898-A-3689.datastack'
tail_clip = 0.05
xlim = [-1, 24]
#def plot_cube_slices(datastack_file, ax=None, tail_clip=0.05, xlim=[-1, 24]):
# read data
ds = DataStack(datastack_file)
x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum')
slices, y_gauss_list = get_cube_slices(datastack_file, tail_clip=tail_clip)
n_slices = len(slices)
fig, ax = plt.subplots(figsize=[9, 4])
# tab20x2 color map
# (there should be an easier way to cycle colors)
tab20 = cm.tab20(np.arange(20))[:,0:3]
tab20_r = cm.tab20_r(np.arange(20))[:,0:3]
tab20x2 = np.r_[tab20, tab20[2:]**0.4]
colors = tab20x2[0:n_slices]
# color gaussian peaks
for i, y_gauss in enumerate(y_gauss_list):
ax.fill_between(x_keVs, y_gauss, color=colors[i])
# color corresponding slice
for i, [si, sj, sk] in enumerate(slices):
is_slice = np.zeros_like(x_keVs)
is_slice[si:sk+1] = 1
y_max = y_gauss_list[i].max() * np.ones_like(x_keVs)
ax.fill_between(x_keVs, y_max, where=is_slice, color=colors[i], alpha=0.3)
ax.plot(x_keVs, y_max, color='r', alpha=0.5, label='max spectrum')
#ax.fill_between(x_keVs, y_max, color='r', alpha=0.2)
_add_hotlines_ticklabels(datastack_file, ax)
ax.set_xlim(xlim)
ax.set_xlabel('energy [keV]')
ax.set_ylabel('intensity [#counts]')
ax.legend()
plt.tight_layout()
# return ax
plot_cube_slices('RP-T-1898-A-3689.datastack')
def plot_cube_slices(datastack_file, ax=None, tail_clip=0.05, xlim=[-1, 24]):
# read data
ds = DataStack(datastack_file)
x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum')
if ax is None:
fig, ax = plt.subplots(figsize=[9, 4])
y_gauss_list = get_cube_slices(datastack_file, tail_clip=tail_clip)[1]
for y_gauss in y_gauss_list:
ax.fill_between(x_keVs, y_gauss)
ax.plot(x_keVs, y_max, color='r', alpha=0.5, label='max spectrum')
ax.fill_between(x_keVs, y_max, color='r', alpha=0.2)
_add_hotlines_ticklabels(datastack_file, ax)
ax.set_xlim(xlim)
ax.set_xlabel('energy [keV]')
ax.set_ylabel('intensity [#counts]')
ax.legend()
plt.tight_layout()
return ax
While trying to extract a smooth scattering baseline from the sum spectrum, my eye was drawn to a peak at the Zinc $K_{\alpha}$ energy. Because we see no copper, it seems that there is a low concentration of zinc present at a large surface in the drawing. I do not understand the source for that? Is is a modern conservation material? In order to find out I need a good peak map. Let's see if we can make one without yet having implemented the baseline estimation and peak overlap corrections...
fig, ax, twax = plot_cube_slices('RP-T-1898-A-3689.datastack')
ax.plot(x_keVs, 10 * y_sum, color=[0.3, 1, 0.3], label='sum spectrum')
ax.legend();
ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies')
si, sj, sk = slices[13] # Zn_Ka slice (left, mid, right) index
# read Zn_Ka slice
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
with ProgressBar():
ZnKa_slice = cube[:,:,si:sk+1].compute() # takes a few seconds...
# integrate
ZnKa_map = ZnKa_slice.sum(axis=2)
# and plot
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])
ax.imshow(ZnKa_map, vmax=150); # need to clip?
ax1.imshow(imvis, extent=extent)
x, y = 880, 555
x2, y2 = 881, 559
ax.scatter(x, y, marker='+', color='r')
ax1.scatter(x, y, marker='+', color='orange')
ax.scatter(x2, y2, marker='+', color='r')
ax1.scatter(x2, y2, marker='+', color='orange')
spectrum_xy = cube[x, y].compute()
spectrum_xy2 = cube[x2, y2].compute()
i, j, k = slices[13] # Zn
Zn_band = np.zeros_like(spectrum_xy)
Zn_band[i:k+1] = 1
fig, ax = plt.subplots(figsize=[9, 4])
ax.fill_between(x_keVs, spectrum_xy, y_compton, where=(spectrum_xy>y_compton), label='xy1')
ax.fill_between(x_keVs, spectrum_xy2, y_compton, where=(spectrum_xy2>y_compton), label='xy2')
#ax.fill_between(x_keVs, spectrum_xy - y_compton, y_compton - y_compton, where=(spectrum_xy>y_compton), label='xy1')
ax.plot(x_keVs, y_compton, color='r')
ax.fill_between(x_keVs, 20 * Zn_band, where=(Zn_band>0), color='g', alpha=0.3)
ax.vlines([x_keVs[j]], 0, 20, color='g')
ax.legend()
_add_hotlines_ticklabels('RP-T-1898-A-3689.datastack', ax)
mos.XFluo('Zn', tube_keV=40).plot(ax=ax);
ax.set_title('Hedgehawk plot for two bright Zn Ka pixels')
plt.tight_layout()
mos.XFluo('Zn', tube_keV=40).plot(ax=ax);
si, sj, sk = slices[21] # Rh_Ka Compton
# read Compton_Rh_Ka slice
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
with ProgressBar():
Compton_slice = cube[:,:,si:sk+1].compute() # takes a few seconds...
# integrate
Compton_map = Compton_slice.sum(axis=2)
fig, [ax, ax1, ax2] = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=[9, 5])
ax.imshow(Compton_map, vmin=650); # need to clip?
ax.set_title('Compton')
ax.scatter(*Compton_max[::-1], color='r')
ax1.imshow(imvis, extent=extent)
ax1.set_title('vis')
ax2.imshow(ZnKa_map, vmin=100, vmax=150); # need to clip?
ax2.set_title('Zinc Ka');
Compton_max = np.argwhere(Compton_map==Compton_map.max()).flatten()
Compton_max